This geospatial statistical model uses routinely collected malaria case data, population data and remotely sensed data, such as open and vegetated water bodies, to estimate population living around open water bodies, expected malaria cases, and standardised morbidity ratio (SMR) of malaria. And ultimately, quantify the association between proximity to larval habitat and malaria risk in health facility catchment areas in Kasungu. The SMR compares the risk of morbidity in a population of interest with that of a standard population. In this case, our interest is to find out whether the number of dry season malaria cases in each catchment area are greater than we would expect given the malaria rate for the entire Kasungu district.

We do this by comparing what we observe (O) with what we would expect (E) if the risk of malaria was equal throughout Kasungu. The SMR statistical notation of catchment i can be written as follows: \[SMR_i = \frac{O_i}{E_i}\]

Buffers around waterbodies are created and then combined with population data in raster format to estimate the proprtion of catcment population living within 1km, 2km and 3km of water bodies. Subsequently, the observed malaria cases are modeled using Poisson regression to find out if living within various distances from water bodies is causing variability in malaria risk in Kasungu district. We hypothesize that the risk of being a case in a catchment is dependent on proximity to water bodies. The data used spans from 2017 to 2020 and was derived from digitized DHIS2 malaria records, accessibility mapping, aggregated population geospatial layer and TropWet tool in Google Earth Engine.

Load packages

Loading the R packages that will be used to read in, view, transform and model the malaria cases and spatial datasets.

library(spatialEco)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(plotly)
library(lubridate)
library(knitr)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(sp)
library(tmap)
library(spdep)
library(maptools)
library(gridExtra)
library(grid)
library(exactextractr)
library(DataExplorer)
library(mapview)
`%>%` <- magrittr::`%>%`

Tell R where the data is

here::here()
## [1] "C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021"

Load datasets

The total dry season malaria cases recorded at health-care facilities in Kasungu from 2017 to 2019 are contained in the KasunguData.csv sourced from https://dhis2.health.gov.mw/. The kasungu_facility_catchments_2004.shp shapefile also contains the population and health information within each health-facility catchment area in Kasungu district.

The aggregated population raster layers for Malawi e.g.,ku_pop_2017_1km_aggregated.tif were downloaded from the Open Spatial and Demographic and Data Research website: https://www.worldpop.org/geodata/country?iso3=MWI. These layers estimate total number of people per grid-cell. The units are number of people per pixel with country totals adjusted to match the corresponding official United Nations population estimates. The datasets were downloaded in Geotiff at a resolution of 1km and are projected in Geographic Coordinate System, WGS84.

The kasungu_water.shpand water_bodies layers contain open and vegetated waterbodies polygons, detected using the Tropical Wetland Unmixing Tool (TropWet). TropWet is a Google Earth Engine hosted toolbox that uses the Landsat archive to map tropical wetlands and can be accessed through: https://www.aber.ac.uk/en/dges/research/earth-observation-laboratory/research/tropwet/

# 2017, 2018 and 2019 dry season malaria cases 
dry_season_malaria_2015_2019 <- read.csv(here::here("data", "dry_season_malaria 2015-2019.csv"))

# Export 2017, 2018 and 2019 dry season malaria cases
write.csv(dry_season_malaria_2015_2019, file = "data/dry_season_malaria_2017_2019.csv")

write.table(dry_season_malaria_2015_2019, 
            file = "dry_season_malaria_2017_2019.txt", 
            sep = ",", 
            quote = FALSE, 
            row.names = FALSE)

# 2020 dry season malaria cases
ku_malaria_2020 <- read.csv(here::here("data", "dry_season_malaria_2020.csv"))


# Merge 2015 to 2019 dry season malaria case data with 2020 data
dry_season_malaria_2017_2020 <- cbind.data.frame(dry_season_malaria_2015_2019, ku_malaria_2020) 
 
dry_season_malaria_2017_2020 <- dry_season_malaria_2017_2020[,c("FID", "Names", "dr_2017", 
                                                                "dr_2018", "dr_2019", "dr_2020",
                                                                "LONGITU", "LATITUD")]

# Export 'dry season malaria 2017-2020' as csv
write.csv(dry_season_malaria_2017_2020, file = "data/dry_season_malaria_2017_2019.csv")


# Kasungu district boundary shapefile 
kasungu_district <- sf::st_read(here::here("data", "kasungu_district.shp"))
## Reading layer `kasungu_district' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_district.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 491272.7 ymin: 8494349 xmax: 609044.2 ymax: 8632164
## Projected CRS: WGS 84 / UTM zone 36S
# Kasungu health facility catchments generated from accessibility mapping

malire_new <- sf::st_read(here::here("data", "new_catchments.shp")) %>% 
              sf::st_transform(32736) # reproject to WGS UTM Zone 36 South
## Reading layer `new_catchments' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\new_catchments.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.925 ymin: -13.61667 xmax: 34.00833 ymax: -12.375
## Geodetic CRS:  WGS 84
# Kasungu population raster layer
kasungu_population_2017 <- raster(here::here("data", "ku_pop_2017_1km_aggregated.tif"))

kasungu_population_2018 <- raster(here::here("data", "ku_pop_2018_1km_aggregated.tif"))

kasungu_population_2019 <- raster(here::here("data", "ku_pop_2019_1km_aggregated.tif"))

kasungu_population_2020 <- raster(here::here("data", "ku_pop_2020_1km_aggregated.tif"))

# Read in waterbodies polygons 
dryseason_waterbodies_2017 <- sf::st_read(here::here("data", "water_bodies_2017.shp"))
## Reading layer `water_bodies_2017' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2017.shp' using driver `ESRI Shapefile'
## Simple feature collection with 168 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 514497 ymin: 8495941 xmax: 603149.8 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36S
dryseason_waterbodies_2018 <- sf::st_read(here::here("data", "kasungu_2018_water.shp"))
## Reading layer `kasungu_2018_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2018_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1105 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 496807.6 ymin: 8494693 xmax: 607913.8 ymax: 8607747
## Projected CRS: WGS 84 / UTM zone 36S
dryseason_waterbodies_2019 <- sf::st_read(here::here("data", "kasungu_2019_water.shp"))
## Reading layer `kasungu_2019_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2019_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1941 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 494197.2 ymin: 8494693 xmax: 607913.8 ymax: 8617573
## Projected CRS: WGS 84 / UTM zone 36S
dryseason_waterbodies_2020 <- sf::st_read(here::here("data", "water_bodies_2020.shp"))
## Reading layer `water_bodies_2020' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 266 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 508985.6 ymin: 8495793 xmax: 585761.1 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36S
# Add a field ID to water bodies polygons 
dryseason_waterbodies_2017$ID <- 1:nrow(dryseason_waterbodies_2017)

dryseason_waterbodies_2018$ID <- 1:nrow(dryseason_waterbodies_2018)

dryseason_waterbodies_2019$ID <- 1:nrow(dryseason_waterbodies_2019)

dryseason_waterbodies_2020$ID <- 1:nrow(dryseason_waterbodies_2020)

View the dry season malaria case data

We observe that Kasungu district has 30 health facilities classified as dispensary, health centre, district hospital and rural hospital, and the highest malaria cases were recorded at Kasungu District Hospital.

dry_season_malaria_2017_2020 %>% 
  summary()
##       FID           Names              dr_2017         dr_2018     
##  Min.   : 1.00   Length:36          Min.   :  427   Min.   :  749  
##  1st Qu.: 9.75   Class :character   1st Qu.: 2196   1st Qu.: 2424  
##  Median :18.50   Mode  :character   Median : 3132   Median : 3357  
##  Mean   :18.50                      Mean   : 3586   Mean   : 3909  
##  3rd Qu.:27.25                      3rd Qu.: 3993   3rd Qu.: 4684  
##  Max.   :36.00                      Max.   :16289   Max.   :15821  
##                                     NA's   :6       NA's   :6      
##     dr_2019         dr_2020         LONGITU         LATITUD      
##  Min.   :  533   Min.   :    0   Min.   :33.18   Min.   :-13.57  
##  1st Qu.: 1748   1st Qu.: 1650   1st Qu.:33.38   1st Qu.:-13.25  
##  Median : 2556   Median : 2698   Median :33.50   Median :-12.98  
##  Mean   : 2880   Mean   : 3404   Mean   :33.52   Mean   :-12.99  
##  3rd Qu.: 3284   3rd Qu.: 4722   3rd Qu.:33.68   3rd Qu.:-12.79  
##  Max.   :10721   Max.   :16435   Max.   :33.87   Max.   :-12.42  
##  NA's   :6                       NA's   :6       NA's   :6
dry_season_malaria_2017_2020 %>%  
  plotly::plot_ly(y = ~Names,
                  x = ~dr_2017,
                  type = "bar",
                  orientation = 'h',
                  name = "2017") %>%
  plotly::add_trace(x = ~ dr_2018,
                    name = "2018") %>%
  plotly::add_trace(x = ~ dr_2019,
                    name = "2019") %>% 
  plotly::add_trace(x = ~ dr_2020,
                    name = "2020") %>% 
  plotly::layout(xaxis = list(title = "Total malaria cases"),
                 yaxis = list(title = " "),
                 hovermode = "compare",
                 margin = list(b = 10,
                               t = 10,
                               pad = 2))

Fig.1 The total malaria cases recorded at each health-care facility in Kasungu district

Kasungu health-care facilities and their catchment areas

Heath facility catchment area is the area from which a health facility attracts patients. The new health facility catchments polygon was generated from generic accessibility mapping script adapted from https://malariaatlas.org/wp-content/uploads/accessibility/R_generic_accessibilty_mapping_script.r The script requires two user supplied datasets: the 2015 friction surface, which is available here: http://www.map.ox.ac.uk/accessibility_to_cities/, and a user-supplied .csv of points dry_season_malaria_2017_2020. The accumulated cost algorithm accCost and r.Cost algorithm in GRASS GIS were run to make the final output map of new health facility catchment boundaries.

# Using the complete.cases() function to select health centres with complete 
# longitude and latitude coordinates.
zipatala_aggregated <- dry_season_malaria_2017_2020[complete.cases(dry_season_malaria_2017_2020),] 

# Aggregate health facilities close to each other: 
#  a) Kasalika Health Centre and Kasungu District Hospital, 
#  b) Bua and Mziza Health Centres, and 
#  c) Kaluluma and Nkhamenya Rural Hospitals in order to 
# generate catchment areas that are geographically correct

zipatala_aggregated$dr_2017[which(
  zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Kasalika Health Centre")]

zipatala_aggregated$dr_2018[which(
  zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2018[which(
    zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2018[which(
      zipatala_aggregated$Names == "Kasalika Health Centre")]

zipatala_aggregated$dr_2019[which(
  zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2019[which(
    zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2019[which(
      zipatala_aggregated$Names == "Kasalika Health Centre")]

zipatala_aggregated$dr_2020[which(
  zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2020[which(
    zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2020[which(
      zipatala_aggregated$Names == "Kasalika Health Centre")]
  
zipatala_aggregated$dr_2017[which(
  zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Kaluluma Rural Hospital")]

zipatala_aggregated$dr_2018[which(
  zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2018[which(
    zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] +zipatala_aggregated$dr_2018[which(
      zipatala_aggregated$Names == "Kaluluma Rural Hospital")]

zipatala_aggregated$dr_2019[which(
  zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2019[which(
    zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2019[which(
      zipatala_aggregated$Names == "Kaluluma Rural Hospital")]

zipatala_aggregated$dr_2020[which(
  zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2020[which(
    zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2020[which(
      zipatala_aggregated$Names == "Kaluluma Rural Hospital")]

zipatala_aggregated$dr_2017[which(
  zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Bua Health Centre")]

zipatala_aggregated$dr_2018[which(
  zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Bua Health Centre")]

zipatala_aggregated$dr_2019[which(
  zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Bua Health Centre")] 
  

zipatala_aggregated$dr_2020[which(
  zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
    zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
      zipatala_aggregated$Names == "Bua Health Centre")]

# Drop out the other health facilities
zipatala_aggregated <- zipatala_aggregated %>%
  dplyr::filter(Names != "Kasalika Health Centre",
                Names != "Bua Health Centre",
                Names != "Kaluluma Rural Hospital")

# write.csv(zipatala_aggregated, "data/health_facilities_aggregated.csv")

# Convert to csv spatial to spatial object
health_facility_aggr_sf <- sf::st_as_sf(zipatala_aggregated,
                                        coords = c("LONGITU", "LATITUD"),
                                        crs = 4326, agr = "constant")

# st_write(health_facility_aggr_sf, "data/health_facilities_aggregated.shp")

View location of the health facilities in the new catchment areas

# Plot map
tm_shape(malire_new)+
  tm_polygons()+
  tm_shape(health_facility_aggr_sf)+
  tm_dots(size = .3, 
          col = "blue", 
          alpha = 0.5)+
  tm_text("Names", 
          size = .3, 
          just = "top", 
          col = "black", 
          remove.overlap = TRUE)+
  tm_layout(frame = FALSE,
            title = "New Kasungu health facility \n catchment boundaries",
            title.size = .8, 
            title.position = c("left", "top"))+
  tm_compass(position=c("right", "top"))+
  tm_scale_bar(breaks = c(0, 10, 20), 
               text.size = .5)
Fig 2. Kasungu health-care facilities and catchment areas

Fig 2. Kasungu health-care facilities and catchment areas

Kasungu district estimated population per grid-cell

# Take a glimpse at the WorldPop raster layers
kasungu_population_2017
## class      : RasterLayer 
## dimensions : 150, 128, 19200  (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667  (x, y)
## extent     : 491272.7, 609044.2, 8494349, 8632164  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2017_1km_aggregated.tif 
## names      : ku_pop_2017_1km_aggregated
kasungu_population_2018
## class      : RasterLayer 
## dimensions : 150, 128, 19200  (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667  (x, y)
## extent     : 491272.7, 609044.2, 8494349, 8632164  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2018_1km_aggregated.tif 
## names      : ku_pop_2018_1km_aggregated 
## values     : 0, 6253.557  (min, max)
kasungu_population_2019
## class      : RasterLayer 
## dimensions : 150, 128, 19200  (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667  (x, y)
## extent     : 491272.7, 609044.2, 8494349, 8632164  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2019_1km_aggregated.tif 
## names      : ku_pop_2019_1km_aggregated 
## values     : 0, 6483.727  (min, max)
kasungu_population_2020
## class      : RasterLayer 
## dimensions : 150, 128, 19200  (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667  (x, y)
## extent     : 491272.7, 609044.2, 8494349, 8632164  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2020_1km_aggregated.tif 
## names      : ku_pop_2020_1km_aggregated 
## values     : 0, 7949.033  (min, max)
# Helper function to create a raster population map
create.population.map <- function(population.raster, title){
  # raster population map
  # arguments:
  #   population.raster:  aggregated population raster layer from WorldPop
  #   legend.title: legend title
  # returns:
  #   a tmap-element (plots a map)
  tm_shape(population.raster)+
    tm_raster(palette = "YlOrBr", 
              title = title,
              breaks = c(0,100,200,400,600,800,1000,2000,4000,6000,8000))+
    tm_layout(legend.position = c("right", "bottom"),
              frame = FALSE)+
    tm_scale_bar(position = c("left", "bottom"))
}
# Set to static map
tmap_mode("plot")

estimated_pop_2017 <- create.population.map(kasungu_population_2017, title = "2017 Population")

estimated_pop_2018 <- create.population.map(kasungu_population_2018, title = "2018 Population")

estimated_pop_2019 <- create.population.map(kasungu_population_2019, title = "2019 Population")

estimated_pop_2020 <- create.population.map(kasungu_population_2020, title = "2020 Population")

# Layout the maps
tmap_arrange(estimated_pop_2017, estimated_pop_2018, estimated_pop_2019, estimated_pop_2020, nrow = 2) 
Fig.3 Estimated total number of people per 1km grid-cell

Fig.3 Estimated total number of people per 1km grid-cell

Assign dry season malaria cases and population density to new health facility catchments

The WorldPop aggregated population e.g. kasungu_population_2017.tif, and DHIS2 malaria dry_season_malaria_2017_2020 datasets are assigned to the new health facility catchments.

# Helper function that assigns malaria data from health facilities to their catchments areas ----------------
assign.malaria.data <- function(catchment_boundary, malaria_data){
  # arguments:
  #   catchment_boundary: sf polygon object of new catchment boundaries
  #   malaria_data: sf point object with a data frame containing the dry season malaria cases
  # returns:
  #   catchments_malaria_sf: sf polygon object with a data frame containing dry season malaria cases


  # Convert sf objects to spatial
  catchment_shp <- as(catchment_boundary, "Spatial")
  
  malaria_shp <- as(malaria_data, "Spatial")

  # Match CRS
  malaria_shp <- spTransform(malaria_shp, crs(catchment_shp))

  # Overlay aggregated health facility points and extract 2017 - 2020 malaria cases
  # Using 'point.in.poly' to return a point spatial object, in this case location of health facilities
  # and estimated population instead of sp::over function, which simply returns 
  # a data frame, with the same no. rows.
  # Argument 'sp = TRUE' returns an sp class object, else returns sf class object
  # Joining the malaria and population dataset using only 'merge' function can't work due to 
  # non-unique columns and differences in row numbers
  
  hospitals_in_catchment <- spatialEco::point.in.poly(malaria_shp, catchment_shp, sp = TRUE) 

  # Add the extracted ID, health facility names and dry season malaria cases to 
  # the health facility catchments (hfc)
  hfc_malaria_shp <- merge(catchment_shp, hospitals_in_catchment, by.x = "DN", by.y = "FID")

  # Convert the shapefile containing malaria data to sf-object
  hfc_malaria_sf <- sf::st_as_sf(hfc_malaria_shp)

  # Tidy the data by dropping columns not needed
  catchment_malaria <- hfc_malaria_sf %>% 
    dplyr::select(-c(coords.x1, coords.x2))

  return(out = catchment_malaria)
}


# Invoking the function ----------------------------------------------------------------------------------
malaria_by_catchment <- assign.malaria.data(malire_new, health_facility_aggr_sf)

Assign population data to the health catchment areas

# Helper unction to extract population from WorldPop raster file and assign ---------------------------
# the values to the new catchments.

extract.pop.values <- function(kasungu_pop_raster, catchments){
  # function to extract population from raster file and assign the population to catchments
  # arguments:
  #   kasungu_pop_raster: population raster file clipped to Kasungu district
  #   catchments: shapefile containing the polygons that we wish to use as boundaries
  # returns:
  #   catchments_malaria_pop_sf: sf polygon object containing malaria and population data
  
  # convert from sf to sp
  catchments_sp <- as(catchments, "Spatial")
  
  # Match extent i.e projection
  catchments_sp <- spTransform(catchments_sp, proj4string(kasungu_pop_raster))
  
  # Crop and mask the population raster to exclude Kasungu National Park
  pop_raster_clip <- raster::mask(raster::crop(kasungu_pop_raster, extent(catchments_sp)), catchments_sp)
  
  # Extracting zonal statistics from a population raster layer. 
  # The population raster is a continuous gridded surface layer that has an 
  # estimated population density value to every square in their grid. 
  # The population values are then summed and apportioned to the catchment polygons
  # catchments_malaria_pop <- catchments %>% 
  #   dplyr::mutate(pop = round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE)))
  
  pop_by_catchment <- round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE))
  
  pop_by_catchment_df <-  pop_by_catchment %>%  
  # apply unlist to the lists to have vectors as the list elements
  lapply(unlist) %>% 
  # convert vectors to data.frames
  lapply(as_tibble) %>% 
  # combine the list of data.frames
  bind_rows(., .id = "rowID") %>% 
  # rename the value variable
  dplyr::rename(pop = value)
  
  # Add row ID to column to catchment layer
  catchments$rowID <- 1:nrow(catchments)
  
  # Merge catchment areas and population data 
  pop_by_catchments <- merge(catchments, pop_by_catchment_df, by.x = "rowID", by.y = "rowID")
  
  # Cleaning 'Inf' values
  pop_by_catchments %>% 
    dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>% 
    dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))

  return(out = pop_by_catchments)
  
}

# Invoking the function ---------------------------------------------------------------------------------------
malaria_pop_by_catchment_2017 <- extract.pop.values(kasungu_population_2017, malaria_by_catchment)

malaria_pop_by_catchment_2018 <- extract.pop.values(kasungu_population_2018, malaria_by_catchment)

malaria_pop_by_catchment_2019 <- extract.pop.values(kasungu_population_2019, malaria_by_catchment)

malaria_pop_by_catchment_2020 <- extract.pop.values(kasungu_population_2020, malaria_by_catchment)

View Kasungu population by catchment maps

Estimated total number of people within health facility catchment areas.

# Helper function to create maps of estimated population by catchment areas --------------------------------
create.population.map <- function(catchment.area, 
                                  variable = "pop", 
                                  title, 
                                  legend.title = "Estimated \n population"){
  # estimated population map
  # catchment.area: estimated population layer from nachulu function
  # variable: variable name (as character, in qoutes)
  # title: map title in quotes
  # legend.title: legend title in qoutes
  # returns:
  #   a tmap-element (plots a map)
  tm_shape(catchment.area)+
    tm_fill(col = variable, 
            breaks = c(0, 13000, 19000, 27000, 35000, 70000, 140000, 200000),
            palette = "YlOrBr",
            title = legend.title)+
    tm_borders(col = "grey",
               lwd = 0.4)+
    tm_layout(legend.position = c(0.75, "bottom"),
              legend.text.size = 0.6,
              legend.title.size = 0.8,
              frame = FALSE)+
    tm_credits(title, 
               position = c(0.3, 0.8), 
               size = 1)
}

# Invoking the function --------------------------------------------------------------------------------
pop_by_catchment_2017 <- create.population.map(malaria_pop_by_catchment_2017, title = "2017")

pop_by_catchment_2018 <- create.population.map(malaria_pop_by_catchment_2018, title = "2018")

pop_by_catchment_2019 <- create.population.map(malaria_pop_by_catchment_2019, title = "2019")

pop_by_catchment_2020 <- create.population.map(malaria_pop_by_catchment_2020, title = "2020")

tmap::tmap_arrange(pop_by_catchment_2017, pop_by_catchment_2018,
                   pop_by_catchment_2019, pop_by_catchment_2020, ncol = 2)
Fig. 4: Estimated population by health facility catchment areas

Fig. 4: Estimated population by health facility catchment areas

Population density by catchment

# Helper function to calculate population density by catchment -----------------
calculate.population.density <- function(pop.data){
  
  # Convert to spatial object
  pop.sp <- as(pop.data, "Spatial")

  # Calculate area of catchment polygon in square kilometres
  pop.sp$area_sqkm <- round(rgeos::gArea(pop.sp, byid = TRUE) / (1000 * 1000))

  # Calculate population density
  pop.sp$pop_density <- round(pop.sp$pop / pop.sp$area_sqkm)
  
  # Convert back to sf object
  pop.sf <- sf::st_as_sf(pop.sp)
  
  return(pop.sf)
}

# Invoking function ------------------------------------------------------------
pop_density_2017 <- calculate.population.density(malaria_pop_by_catchment_2017)

pop_density_2018 <- calculate.population.density(malaria_pop_by_catchment_2018)

pop_density_2019 <- calculate.population.density(malaria_pop_by_catchment_2019)

pop_density_2020 <- calculate.population.density(malaria_pop_by_catchment_2020)

# Helper function to create population density maps ----------------------------
create.pop.density.map <- function(pop.density.data,
                                   variable = "pop_density", 
                                   title = NA, 
                                   legend.title = "Population \ndensity/sqkm"){
  tm_shape(pop.density.data)+
    tm_fill(col = variable, 
            breaks = c(0, 50, 100, 150, 200, 250, 300, 350),
            palette = "-magma",
            title = legend.title)+
    tm_borders()+
    tm_layout(legend.position = c(0.75, "bottom"),
              legend.text.size = 0.6,
              legend.title.size = 0.8,
              frame = FALSE)+
    tm_credits(title, 
               position = c(0.3, 0.8), 
               size = 1)
}

# Invoking function ------------------------------------------------------------
pop_density_2017_map <- create.pop.density.map(pop_density_2017, title = "2017")

pop_density_2018_map <- create.pop.density.map(pop_density_2018, title = "2018")

pop_density_2019_map <- create.pop.density.map(pop_density_2019, title = "2019")

pop_density_2020_map <- create.pop.density.map(pop_density_2020, title = "2020")

# Layout maps
tmap::tmap_arrange(pop_density_2017_map, pop_density_2018_map,
                   pop_density_2019_map, pop_density_2020_map, ncol = 2)
Fig. 5: Estimated population density by health facility catchment areas

Fig. 5: Estimated population density by health facility catchment areas

Calculate the expected number of cases for each catchment area

The expected number of dry season malaria cases in catchment i are calculated as the observed risk (r) of malaria i.e. the total number of malaria cases in Kasungu district divided by the total population of the district, multiplied by the number of people in the catchment area: \[E_i = \frac{\sum_i O_i}{\sum_i N_i}\times N_i\]

The expected number of dry season malaria cases are calculated under the assumption that there is no spatial variation in risk, i.e., no difference in infection rates between the catchment areas.

# Calculate expected malaria cases --------------------------------------------------------------
expected_malaria_2017 <- malaria_pop_by_catchment_2017 %>% 
  dplyr::rename(
    observed_2017 = dr_2017,
     pop_2017 = pop) %>% 
  dplyr::mutate(
    expected_2017 = round(sum(observed_2017)/sum(pop_2017, na.rm = TRUE)*pop_2017))

expected_malaria_2018 <- malaria_pop_by_catchment_2018 %>% 
  dplyr::rename(
    observed_2018 = dr_2018,
    pop_2018 = pop) %>% 
  dplyr::mutate(
    expected_2018 = round(sum(observed_2018)/sum(pop_2018, na.rm = TRUE)*pop_2018))

expected_malaria_2019 <- malaria_pop_by_catchment_2019 %>% 
  dplyr::rename(
    observed_2019 = dr_2019,
    pop_2019 = pop) %>%
  dplyr::mutate(
    expected_2019 = round(sum(observed_2019)/sum(pop_2019, na.rm = TRUE)*pop_2019)) 

expected_malaria_2020 <- malaria_pop_by_catchment_2020 %>% 
  dplyr::rename(
    observed_2020 = dr_2020,
    pop_2020 = pop) %>% 
  dplyr::mutate(
    expected_2020 = round(sum(observed_2020)/sum(pop_2020, na.rm = TRUE)*pop_2020))

Calculate the Standardised Morbidity Ratio of malaria incidences for each catchment area

The SMR compares the risk of morbidity in a population of interest with that of a standard population. In this case, our interest is to find out whether the number of dry season malaria cases in each catchment area are greater than we would expect given the malaria rate for the entire Kasungu district.

We do this by comparing what we observe (O) with what we would expect (E) if the risk of malaria was equal throughout Kasungu. The SMR statistical notation of catchment i can be written as follows: \[SMR_i = \frac{O_i}{E_i}\]

# Calculate Standardised Morbidity Ratio (SMR) -------------------------------------
SMR_2017 <- expected_malaria_2017 %>% 
  dplyr::mutate(SMR = round(observed_2017/expected_2017, 1)) %>% 
  dplyr::select(rowID,Names, pop_2017, observed_2017, expected_2017, SMR) 

SMR_2018 <- expected_malaria_2018 %>% 
  dplyr::mutate(SMR = round(observed_2018/expected_2018, 1)) %>% 
  dplyr::select(rowID, Names, pop_2018, observed_2018, expected_2018, SMR) 

SMR_2019 <- expected_malaria_2019 %>% 
  dplyr::mutate(SMR = round(observed_2019/expected_2019, 1)) %>% 
  dplyr::select(rowID, Names, pop_2019, observed_2019, expected_2019, SMR) 

SMR_2020 <- expected_malaria_2020 %>% 
  dplyr::mutate(SMR = round(observed_2020/expected_2020, 1)) %>% 
  dplyr::select(rowID, Names, pop_2020, observed_2020, expected_2020, SMR)


# Create SMR tables ------------------------------------------------------------
SMR_table_2017 <- SMR_2017 %>% 
  dplyr::as_tibble() %>% 
  dplyr::select(-rowID, -geometry) %>% 
  kable %>%
  kableExtra::kable_styling(full_width = FALSE)

SMR_table_2018 <- SMR_2018 %>% 
  dplyr::as_tibble() %>% 
  dplyr::select(-rowID, -geometry) %>% 
  kable %>% 
  kableExtra::kable_styling(full_width = FALSE)

SMR_table_2019 <- SMR_2019 %>% 
  dplyr::as_tibble() %>% 
  dplyr::select(-rowID, -geometry) %>% 
  kable %>% 
  kableExtra::kable_styling(full_width = FALSE)

SMR_table_2020 <- SMR_2020 %>% 
  dplyr::as_tibble() %>% 
  dplyr::select(-rowID, -geometry) %>% 
  kable %>% 
  kableExtra::kable_styling(full_width = FALSE)

SMR_table_2017
Names pop_2017 observed_2017 expected_2017 SMR
Lodjwa Health Centre 9923 1161 1372 0.8
Nkhamenya Rural Hospital 40154 3639 5553 0.7
Newa Mpasazi Health Centre 13879 427 1919 0.2
Mpepa /Chisinga Health Centre 27459 2784 3797 0.7
Mnyanja Health Centre 39950 3298 5524 0.6
Simlemba Health Centre 26999 2197 3734 0.6
Ofesi Health Centre 28098 4036 3886 1.0
Chulu Health Centre 27906 5801 3859 1.5
Kapelula Health Centre 35727 4138 4940 0.8
Livwezi Health Centre 22009 1307 3043 0.4
Gogode Dispensary 13061 2745 1806 1.5
Dwangwa Dispensary 32704 2192 4522 0.5
Chamama Health Facility 20026 1878 2769 0.7
Wimbe Health Centre 11864 5660 1641 3.4
Chinyama 12768 2292 1766 1.3
Mdunga Health Centre 18177 2045 2514 0.8
Mtunthama Health Centre 18744 3622 2592 1.4
Kasungu District Hospital 143490 20817 19842 1.0
Chamwabvi Health Centre 35353 3601 4889 0.7
Linyangwa Health Centre 17772 3359 2458 1.4
Kawamba Health Centre 22865 5808 3162 1.8
Mziza Health Centre 44189 7352 6111 1.2
Kamboni Health Centre 21226 3824 2935 1.3
Khola Health Centre 16956 2195 2345 0.9
Santhe Health Centre 6096 5838 843 6.9
Anchor Farm 48861 2966 6757 0.4
Mkhota Health Centre 21621 2586 2990 0.9
SMR_table_2018
Names pop_2018 observed_2018 expected_2018 SMR
Lodjwa Health Centre 10281 1151 1588 0.7
Nkhamenya Rural Hospital 41642 4540 6430 0.7
Newa Mpasazi Health Centre 14248 749 2200 0.3
Mpepa /Chisinga Health Centre 28488 3602 4399 0.8
Mnyanja Health Centre 41856 2864 6464 0.4
Simlemba Health Centre 27455 2515 4240 0.6
Ofesi Health Centre 29002 3446 4479 0.8
Chulu Health Centre 28832 4502 4452 1.0
Kapelula Health Centre 37630 5338 5811 0.9
Livwezi Health Centre 22544 2361 3481 0.7
Gogode Dispensary 13368 4138 2064 2.0
Dwangwa Dispensary 33534 2394 5178 0.5
Chamama Health Facility 20372 1750 3146 0.6
Wimbe Health Centre 11814 5010 1824 2.7
Chinyama 13138 2116 2029 1.0
Mdunga Health Centre 18928 2923 2923 1.0
Mtunthama Health Centre 19074 5308 2945 1.8
Kasungu District Hospital 147175 20314 22727 0.9
Chamwabvi Health Centre 36167 4027 5585 0.7
Linyangwa Health Centre 18032 3268 2785 1.2
Kawamba Health Centre 22902 6462 3537 1.8
Mziza Health Centre 46208 10841 7136 1.5
Kamboni Health Centre 21430 4745 3309 1.4
Khola Health Centre 17315 2802 2674 1.0
Santhe Health Centre 6244 7267 964 7.5
Anchor Farm 49871 3142 7701 0.4
Mkhota Health Centre 22167 5920 3423 1.7
SMR_table_2019
Names pop_2019 observed_2019 expected_2019 SMR
Lodjwa Health Centre 10608 1713 1200 1.4
Nkhamenya Rural Hospital 43293 4537 4896 0.9
Newa Mpasazi Health Centre 14780 809 1672 0.5
Mpepa /Chisinga Health Centre 29456 4248 3331 1.3
Mnyanja Health Centre 43783 3148 4952 0.6
Simlemba Health Centre 28076 2339 3175 0.7
Ofesi Health Centre 30065 4079 3400 1.2
Chulu Health Centre 29731 5470 3362 1.6
Kapelula Health Centre 39747 2558 4495 0.6
Livwezi Health Centre 22945 787 2595 0.3
Gogode Dispensary 13641 2200 1543 1.4
Dwangwa Dispensary 34415 2553 3892 0.7
Chamama Health Facility 20701 1508 2341 0.6
Wimbe Health Centre 11855 3041 1341 2.3
Chinyama 13475 1770 1524 1.2
Mdunga Health Centre 19960 1008 2257 0.4
Mtunthama Health Centre 19385 2763 2192 1.3
Kasungu District Hospital 151079 13450 17086 0.8
Chamwabvi Health Centre 36899 1686 4173 0.4
Linyangwa Health Centre 18279 3330 2067 1.6
Kawamba Health Centre 23041 4402 2606 1.7
Mziza Health Centre 48340 10841 5467 2.0
Kamboni Health Centre 21509 2770 2433 1.1
Khola Health Centre 17761 2708 2009 1.3
Santhe Health Centre 6435 5210 728 7.2
Anchor Farm 50995 1978 5767 0.3
Mkhota Health Centre 22677 2162 2565 0.8
SMR_table_2020
Names pop_2020 observed_2020 expected_2020 SMR
Lodjwa Health Centre 13081 1788 1603 1.1
Nkhamenya Rural Hospital 53692 8539 6581 1.3
Newa Mpasazi Health Centre 18311 2182 2244 1.0
Mpepa /Chisinga Health Centre 36317 5186 4451 1.2
Mnyanja Health Centre 54649 6117 6698 0.9
Simlemba Health Centre 34240 5310 4197 1.3
Ofesi Health Centre 37240 2323 4564 0.5
Chulu Health Centre 36638 7160 4491 1.6
Kapelula Health Centre 50214 7297 6154 1.2
Livwezi Health Centre 27786 1028 3406 0.3
Gogode Dispensary 16681 2767 2045 1.4
Dwangwa Dispensary 42282 2869 5182 0.6
Chamama Health Facility 25248 635 3095 0.2
Wimbe Health Centre 14367 2233 1761 1.3
Chinyama 16463 1605 2018 0.8
Mdunga Health Centre 25108 3169 3077 1.0
Mtunthama Health Centre 23501 1882 2880 0.7
Kasungu District Hospital 185282 19393 22709 0.9
Chamwabvi Health Centre 45106 1128 5528 0.2
Linyangwa Health Centre 22144 4380 2714 1.6
Kawamba Health Centre 27961 7073 3427 2.1
Mziza Health Centre 60510 10841 7416 1.5
Kamboni Health Centre 25750 4665 3156 1.5
Khola Health Centre 21929 3426 2688 1.3
Santhe Health Centre 7917 4891 970 5.0
Anchor Farm 62633 1665 7677 0.2
Mkhota Health Centre 27830 4592 3411 1.3

View observed and expected dry season malaria cases

# Helper function to create maps of observed and expected dry season malaria cases
create.malaria.map <- function(malaria.data, 
                               variable = NA, 
                               title = NA, 
                               legend.title = NA){
  # observed and expected malaria incidence map
  # malaria.data: data frame containing observed and expected malaria cases
  # variable: variable name (as character, in quotes e.g. "observed")
  # title: map title in quotes
  # legend.title: legend title in quotes
  # returns:
  #   a tmap-element (plots a map)
  tm_shape(malaria.data)+
    tm_fill(col = variable, 
            breaks = c(0, 500, 1000, 2500, 5000, 10000, 15000, 20000, 25000),
            palette = "YlOrRd",
            title = legend.title)+
    tm_borders(lw = 0.3)+
    tm_layout(legend.position = c(0.75,"bottom"),
              legend.text.size = 0.5,
              legend.title.size = 0.7,
              frame = FALSE)+
    tm_credits(title, 
               position = c(0.2, 0.8), 
               size = 1)
}

# Invoking the function
# 2017 observed and expected malaria cases -------------------------------------
observed_malaria_2017_map <- create.malaria.map(malaria_pop_by_catchment_2017, 
                                                variable = "dr_2017",
                                                title = "2017",
                                                legend.title = "Observed malaria")

expected_malaria_2017_map <- create.malaria.map(expected_malaria_2017,
                                                variable = "expected_2017",
                                                title = "2017",
                                                legend.title = "Expected malaria")

# 2018 observed and expected malaria cases -------------------------------------
observed_malaria_2018_map <- create.malaria.map(malaria_pop_by_catchment_2018,
                                                variable = "dr_2018",
                                                title = "2018",
                                                legend.title = "Observed malaria")

expected_malaria_2018_map <- create.malaria.map(expected_malaria_2018,
                                                variable = "expected_2018",
                                                title = "2018",
                                                legend.title = "Expected malaria")

# 2019 observed and expected malaria cases -------------------------------------
observed_malaria_2019_map <- create.malaria.map(malaria_pop_by_catchment_2019,
                                                variable = "dr_2019",
                                                title = "2019",
                                                legend.title = "Observed malaria")

expected_malaria_2019_map <- create.malaria.map(expected_malaria_2019,
                                                variable = "expected_2019",
                                                title = "2019",
                                                legend.title = "Expected malaria")

# 2020 observed and expected malaria cases -------------------------------------
observed_malaria_2020_map <- create.malaria.map(malaria_pop_by_catchment_2020,
                                                variable = "dr_2020",
                                                title = "2020",
                                                legend.title = "Observed malaria")

expected_malaria_2020_map <- create.malaria.map(expected_malaria_2020,
                                                variable = "expected_2020",
                                                title = "2020",
                                                legend.title = "Expected malaria")

# Layout maps ------------------------------------------------------------------
tmap::tmap_arrange(observed_malaria_2017_map, expected_malaria_2017_map,
                   observed_malaria_2018_map, expected_malaria_2018_map, 
                   observed_malaria_2019_map, expected_malaria_2019_map,
                   observed_malaria_2020_map, expected_malaria_2020_map, ncol = 2)
Fig 6: Observed and expected malaria incidence by health facility catchment area, Kasungu

Fig 6: Observed and expected malaria incidence by health facility catchment area, Kasungu

View SMR by catchment

A ratio greater than 1.0 indicates that more malaria cases have occurred than would have been expected, while a ratio less than 1.0 indicates that less cases have occurred.

#   max(SMR_2017$SMR)
# [1] 6.93
# > max(SMR_2018$SMR)
# [1] 7.93
# > max(SMR_2019$SMR)
# [1] 7.71
# > max(SMR_2020$SMR)
# [1] 5.68

# Helper function to create maps of SMR by catchment ----------------------------------
create.smr.map <- function(smr.data, 
                           variable = "SMR", 
                           title = NA, 
                           legend.title = "SMR"){
  # SMR by catchment map
  # smr.data: sf polygon object containing SMR by catchment data
  # variable: variable name (as character, in qoutes)
  # title: map title in quotes
  # legend.title: legend title in qoutes
  # returns:
  #   a tmap-element (plots a map)
 
  tm_shape(smr.data)+
    tm_fill(col = variable, 
            breaks = c(0, 0.5, 1, 1.5, 2, 2.5, 5, 8),
            palette = "-magma",
            title = legend.title)+
    tm_borders(lw = 0.3)+
    tm_layout(legend.position = c(0.75,"bottom"),
              legend.text.size = 0.5,
              legend.title.size = 0.7,
              frame = FALSE)+
    tm_credits(title, 
               position = c(0.2, 0.8), 
               size = 1)
}

# Invoking function -------------------------------------------------------------------
SMR_2017_map <- create.smr.map(SMR_2017, title = "2017")

SMR_2018_map <- create.smr.map(SMR_2018, title = "2018")

SMR_2019_map <- create.smr.map(SMR_2019, title = "2019")

SMR_2020_map <- create.smr.map(SMR_2020, title = "2020")

# Layout maps -------------------------------------------------------------------------
tmap::tmap_arrange(SMR_2017_map, SMR_2018_map, SMR_2019_map, SMR_2020_map, ncol = 2)
Fig. 7: Standardised morbidity ratio of malaria by health facility catchment

Fig. 7: Standardised morbidity ratio of malaria by health facility catchment

Calculate the proportion of the catchment population living within 1km, 2km, 3km of water bodies

First, using st_buffer, we compute 1km, 2km and 3km buffers around dry season water bodies obtained from LandSat satellite imagery using TropWet tool in Google Earth Engine. Then geometry of the buffer features are then combined resulting in resolved internal boundaries to enable extracting population values from WorldPop raster. Finally, we calculate the proportion of people in each catchment area living within water bodies.

# Combine and transform TropWet derived waterbody polygons -------------------------------
surface_waterbodies_2017 <- sf::st_as_sf(
  st_cast(
    st_union(
      st_buffer(dryseason_waterbodies_2017, dist = 30)), "POLYGON"))

surface_waterbodies_2018 <- sf::st_as_sf(
  st_cast(
    st_union(
      st_buffer(dryseason_waterbodies_2018, dist = 30)), "POLYGON"))

surface_waterbodies_2019 <- sf::st_as_sf(
  st_cast(
    st_union(
      st_buffer(dryseason_waterbodies_2019, dist = 30)), "POLYGON"))

surface_waterbodies_2020 <- sf::st_as_sf(
  st_cast(
    st_union(
      st_buffer(dryseason_waterbodies_2020, dist = 30)), "POLYGON"))

# Helper function to compute 1km, 2km and 3km buffers around the water bodies ---------------------

create.waterbody.buffer <- function(waterbody, distance, catchment){
  # function for creating buffers around waterbodies
  # arguments:
  #   waterbody:  waterbody shapefile
  #   distance: buffer distance in meters
  #   catchment: catchment area shapefile
  # returns:
  #   buffered waterbodies 
  
  # Create buffers around water bodies
  buffer_radius <- sf::st_buffer(waterbody, distance)
  
  # Dissolve the buffers
  buffer_union <- sf::st_as_sf(st_cast(st_union(buffer_radius),"MULTIPOLYGON"))
  
  # Assign attributes of the 'catchment' to each of the water bodies. 
   buffer_intersect <- sf::st_intersection(buffer_union, catchment)
  
   buffer_intersect_sf <- sf::st_as_sf(buffer_intersect)
   
  # Convert the MULTIPOLYGON object into several POLYGON objects
   buffer_intersect_polygons <- sf::st_cast(
     sf::st_buffer(buffer_intersect_sf,0.0), "MULTIPOLYGON") %>% 
     sf::st_cast("POLYGON")
  
  # Polygons being seen to be in multiple catchments
   sf::st_intersects(buffer_intersect_polygons, catchment)
  
  # Make the assumption that the attribute is constant throughout the geometry
   sf::st_agr(buffer_intersect_polygons) = "constant"
   
   sf::st_agr(catchment) = "constant"
  
  return(out = buffer_intersect_polygons)
}


# Invoking function
# For 2017 TropWet surface water polygons --------------------------------------------------------
buffer_1km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017, 
                                           distance = 1000, 
                                           catchment = malire_new)

buffer_2km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017, 
                                           distance = 2000, 
                                           catchment = malire_new)

buffer_3km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017, 
                                           distance = 3000,
                                           catchment = malire_new)

# For 2018 TropWet surface water polygons --------------------------------------------------------
buffer_1km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018, 
                                           distance = 1000, 
                                           catchment = malire_new)

buffer_2km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018, 
                                           distance = 2000, 
                                           catchment = malire_new)

buffer_3km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018, 
                                           distance = 3000, 
                                           catchment = malire_new)
 
# For 2019 TropWet surface water polygons ------------------------------------------------------
buffer_1km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019, 
                                           distance = 1000, 
                                           catchment = malire_new)

buffer_2km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019, 
                                           distance = 2000, 
                                           catchment = malire_new)

buffer_3km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019, 
                                           distance = 3000, 
                                           catchment = malire_new)

# For 2020 TropWet surface water polygons ------------------------------------------------------
buffer_1km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020, 
                                           distance = 1000, 
                                           catchment = malire_new)

buffer_2km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020, 
                                           distance = 2000, 
                                           catchment = malire_new)

buffer_3km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020, 
                                           distance = 3000, 
                                           catchment = malire_new)

View the created waterbody buffers

# Map the buffers
create.buffer.map <- function(buffers, boundary = malire_new, title = NA){
  # function for creating buffer map in ggplot
  # arguments:
  #   buffer:  waterbodies buffer polygon layer
  #   boundary: health facility catchment polygons
  #   title: main title
  # returns:
  #   a map-element (plots a map)
  ggplot(data = buffers)+
     geom_sf()+
     geom_sf(data = boundary, 
             fill = NA)+
     theme_void()+
     labs(title = title)
}

# Invoking the function
# For 2017 -------------------------------------------------------------------------------
buffer_1km_2017_map <- create.buffer.map(buffer_1km_2017, title = "2017: 1km Buffers")

buffer_2km_2017_map <- create.buffer.map(buffer_2km_2017, title = "2017: 2km Buffers")

buffer_3km_2017_map <- create.buffer.map(buffer_3km_2017, title = "2017: 3km Buffers")

# For 2018 --------------------------------------------------------------------------------
buffer_1km_2018_map <- create.buffer.map(buffer_1km_2018, title = "2018: 1km Buffers")

buffer_2km_2018_map <- create.buffer.map(buffer_2km_2018, title = "2018: 2km Buffers")

buffer_3km_2018_map <- create.buffer.map(buffer_3km_2018, title = "2018: 3km Buffers")

# For 2019 ---------------------------------------------------------------------------------
buffer_1km_2019_map <- create.buffer.map(buffer_1km_2019, title = "2019: 1km Buffers")

buffer_2km_2019_map <- create.buffer.map(buffer_2km_2019, title = "2019: 2km Buffers")

buffer_3km_2019_map <- create.buffer.map(buffer_3km_2019, title = "2019: 3km Buffers")

# For 2020 --------------------------------------------------------------------------------
buffer_1km_2020_map <- create.buffer.map(buffer_1km_2020, title = "2020: 1km Buffers")

buffer_2km_2020_map <- create.buffer.map(buffer_2km_2020, title = "2020: 2km Buffers")

buffer_3km_2020_map <- create.buffer.map(buffer_3km_2020, title = "2020: 3km Buffers")
 
grid.arrange(buffer_1km_2017_map, buffer_1km_2018_map, buffer_1km_2019_map, buffer_1km_2020_map,
             buffer_2km_2017_map, buffer_2km_2018_map, buffer_2km_2019_map, buffer_2km_2020_map, 
             buffer_3km_2017_map, buffer_3km_2018_map, buffer_3km_2019_map, buffer_3km_2020_map, ncol = 4)
Fig 8. Buffers around dry season waterbodies in Kasungu

Fig 8. Buffers around dry season waterbodies in Kasungu

Extract the population living within waterbody buffers by catchment area

# Helper function to calculate estimated number of people living within waterbody buffers
# in each catchment area
estimate.buffer.pop <- function(catchment.population, buffers, catchment.area){
  
  # Extract population estimates from WorldPop raster
  buffers$buffer_pop <- raster::extract(catchment.population,
                                        buffers, 
                                        fun = sum, 
                                        na.rm = TRUE)
                                               
                                              
  # Find which catchment each polygon belongs to using its centroid - a point dataset 
  # representing the geographic center-points of the polygons 
  buffer_by_catchment <- st_intersection(st_centroid(buffers), catchment.area)
  
  # Notice that the buffer_catchment is comprised of separate POLYGONS (buffer_by_catchment$x). 
  # The first step is to “dissolve” away these POLYGONS into one MULTIPOLYGON. 
  # There is no sf equivalent to the QGIS or ArcMap “dissolve” operation. 
  # Instead we use a combination of group_by and summarize from the dplyr package. 
  # Stats::aggregate from sf package, and dplyr::summarize both do essentially the same.
   buffer_pop_aggregated <- buffer_by_catchment %>% 
     dplyr::group_by(DN) %>%
     dplyr::summarize(
       buffer_pop_aggregated = round(sum(buffer_pop, na.rm = TRUE)))
   
  buffer_pop <- merge(
    catchment.area, st_drop_geometry(
      buffer_pop_aggregated), by = 'DN', all.x = TRUE)
  
  return(out = buffer_pop)
  
}

# Invoking the function and calculating proportion of 
# catchment population living within buffers
# 2017 buffer population -------------------------------------------------------
buffer_pop_1km_2017 <- estimate.buffer.pop(
  kasungu_population_2017, 
  buffer_1km_2017, 
   malaria_pop_by_catchment_2017) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100)) 

buffer_pop_2km_2017 <- estimate.buffer.pop(
  kasungu_population_2017,
  buffer_2km_2017,
  malaria_pop_by_catchment_2017) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100)) 

buffer_pop_3km_2017 <- estimate.buffer.pop(
  kasungu_population_2017,
  buffer_3km_2017,
  malaria_pop_by_catchment_2017) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

# 2018 buffer population -------------------------------------------------------
buffer_pop_1km_2018 <- estimate.buffer.pop(
  kasungu_population_2018,
  buffer_1km_2018,
  malaria_pop_by_catchment_2018) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_2km_2018 <- estimate.buffer.pop(
  kasungu_population_2018,
  buffer_2km_2018,
  malaria_pop_by_catchment_2018) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_3km_2018 <- estimate.buffer.pop(
  kasungu_population_2018,
  buffer_3km_2018,
  malaria_pop_by_catchment_2018) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

# 2019 buffer population -------------------------------------------------------
buffer_pop_1km_2019 <- estimate.buffer.pop(
  kasungu_population_2019,
  buffer_1km_2019,
  malaria_pop_by_catchment_2019) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_2km_2019 <- estimate.buffer.pop(
  kasungu_population_2019,
  buffer_2km_2019,
  malaria_pop_by_catchment_2019) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_3km_2019 <- estimate.buffer.pop(
  kasungu_population_2019,
  buffer_3km_2019,
  malaria_pop_by_catchment_2019) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

# 2020 buffer population -------------------------------------------------------
buffer_pop_1km_2020 <- estimate.buffer.pop(
  kasungu_population_2020,
  buffer_1km_2020,
  malaria_pop_by_catchment_2020) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_2km_2020 <- estimate.buffer.pop(
  kasungu_population_2020,
  buffer_2km_2020,
  malaria_pop_by_catchment_2020) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

buffer_pop_3km_2020 <- estimate.buffer.pop(
  kasungu_population_2020,
  buffer_3km_2020,
  malaria_pop_by_catchment_2020) %>% 
  dplyr::rename(catchment_pop = pop,
                buffer_pop = buffer_pop_aggregated) %>% 
  dplyr::mutate(
    prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))

Mapping proportion of catchment population living within waterbodies

# Helper function to create maps of proportion of people living in proximity ----------
# to water bodies in each catchment area
create.pop.proportion.map <- function(pop.data, 
                                      variable = "prop_buffer_catchment_pop", 
                                      title = NA, 
                                      legend.title = NA){
 
  # pop.data: sf polygon object containing proportion of catchment population 
  #           living within water bodies
  # variable: variable name (as character, in qoutes)
  # title: map title in quotes
  # legend.title: legend title in qoutes
  # returns:
  #   a tmap-element (plots a map)
 
  tm_shape(pop.data)+
    tm_fill(col = variable, 
            breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100),
            palette = "YlOrBr",
            title = legend.title)+
    tm_borders(lw = 0.3)+
    tm_layout(legend.position = c(0.8,"bottom"),
              legend.text.size = 0.5,
              legend.title.size = 0.7,
              frame = FALSE)+
    tm_credits(title, 
               position = c(0.25, 0.75), 
               size = 1)
}

# Invoking function 
# 2017 population proportion ---------------------------------------------------
pop_proportion_1km_2017_map <- create.pop.proportion.map(
  buffer_pop_1km_2017, 
  title = "2017",
  legend.title = "Population within \n1km buffers (%)")

pop_proportion_2km_2017_map <- create.pop.proportion.map(
  buffer_pop_2km_2017, 
  title = "2017",
  legend.title = "Population within \n2km buffers (%)")

pop_proportion_3km_2017_map <- create.pop.proportion.map(
  buffer_pop_3km_2017,
  title = "2017",
  legend.title = "Population within \n3km buffers (%)")

# 2018 population proportion ---------------------------------------------------
pop_proportion_1km_2018_map <- create.pop.proportion.map(
  buffer_pop_1km_2018,
  title = "2018",
  legend.title = "Population within \n1km buffers (%)")

pop_proportion_2km_2018_map <- create.pop.proportion.map(
  buffer_pop_2km_2018,
  title = "2018",
  legend.title = "Population within \n2km buffers (%)")

pop_proportion_3km_2018_map <- create.pop.proportion.map(
  buffer_pop_3km_2018,
  title = "2018",
  legend.title = "Population within \n3km buffers (%)")

# 2019 population proportion ---------------------------------------------------
pop_proportion_1km_2019_map <- create.pop.proportion.map(
  buffer_pop_1km_2019,
  title = "2019",
  legend.title = "Population within \n1km buffers (%)")

pop_proportion_2km_2019_map <- create.pop.proportion.map(
  buffer_pop_2km_2019,
  title = "2019",
  legend.title = "Population within \n2km buffers (%)")

pop_proportion_3km_2019_map <- create.pop.proportion.map(
  buffer_pop_3km_2019,
  title = "2019",
  legend.title = "Population within \n3km buffers (%)")

# 2020 population proportion ---------------------------------------------------
pop_proportion_1km_2020_map <- create.pop.proportion.map(
  buffer_pop_1km_2020,
  title = "2020",
  legend.title = "Population within \n1km buffers (%)")

pop_proportion_2km_2020_map <- create.pop.proportion.map(
  buffer_pop_2km_2020,
  title = "2020",
  legend.title = "Population within \n2km buffers (%)")

pop_proportion_3km_2020_map <- create.pop.proportion.map(
  buffer_pop_3km_2020,
  title = "2020",
  legend.title = "Population within \n3km buffers (%)")

# Layout maps ------------------------------------------------------------------
tmap::tmap_arrange(pop_proportion_1km_2017_map, pop_proportion_2km_2017_map, 
                   pop_proportion_3km_2017_map, pop_proportion_1km_2018_map,
                   pop_proportion_2km_2018_map, pop_proportion_3km_2018_map,
                   pop_proportion_1km_2019_map, pop_proportion_2km_2019_map,
                   pop_proportion_3km_2019_map, pop_proportion_1km_2020_map,
                   pop_proportion_2km_2020_map, pop_proportion_3km_2020_map, ncol = 3)
Fig 8. Proportion of catchment population living around water bodies

Fig 8. Proportion of catchment population living around water bodies

Scatter plots of SMR against the proportion of the catchment population living waterbody buffers

A correlation coeeficient of more than zero (cor.coeff r > 0.1) indicates some positive association between the SMR and the buffer population variables. That is, SMR of dry season malaria increases with increase in number of people suurounding water bodies.

# Helper function to tidy and bind the SMR and proportion of -------------------
# buffer-catchment population data frames
tidy.data <- function(smr.df, 
                      proportion.pop.1km, 
                      proprotion.pop.2km,
                      proportion.pop.3km){

# Convert the sf object to dataframes-------------------------------------------
smr_df <- as.data.frame(smr.df) %>% 
  dplyr::select(rowID, Names, SMR)

proportion_pop_1km_df <- as.data.frame(proportion.pop.1km) %>% 
  dplyr::select(rowID, prop_pop_1km = `prop_buffer_catchment_pop`)

proportion_pop_2km_df <- as.data.frame(proprotion.pop.2km)%>% 
  dplyr::select(rowID, prop_pop_2km = `prop_buffer_catchment_pop`)

proportion_pop_3km_df <- as.data.frame(proportion.pop.3km)%>% 
  dplyr::select(rowID, prop_pop_3km = `prop_buffer_catchment_pop`)

# Merge SMR and population data frames -----------------------------------------
combined_1 <- merge(smr_df, proportion_pop_1km_df, by = "rowID", all = TRUE)

combined_2 <- merge(proportion_pop_2km_df, proportion_pop_3km_df)

combined_fully <- merge(combined_1, combined_2, by = "rowID", all = TRUE)

}

# Invoking the function --------------------------------------------------------
smr_pop_2017 <- tidy.data(SMR_2017, buffer_pop_1km_2017, buffer_pop_2km_2017, buffer_pop_3km_2017)

smr_pop_2018 <- tidy.data(SMR_2018, buffer_pop_1km_2018, buffer_pop_2km_2018, buffer_pop_3km_2018)

smr_pop_2019 <- tidy.data(SMR_2019, buffer_pop_1km_2019, buffer_pop_2km_2019, buffer_pop_3km_2019)

smr_pop_2020 <- tidy.data(SMR_2020, buffer_pop_1km_2020, buffer_pop_2km_2020, buffer_pop_3km_2020)

# Helper function to create scatter plots --------------------------------------
create.scatter.plot <- function(smr.pop.df, 
                                independent.var = NA,
                                dependent.var = "SMR",
                                x.label = NA,
                                plot.title = NA){
  
  scatter.plot <- ggpubr::ggscatter(smr.pop.df,          # data frame
                                    x = independent.var, # x-axis variable
                                    y = dependent.var,   # y-axis variable
                                    add = "reg.line",    # Add regression line
                                    conf.int = TRUE,     # Add confidence interval
                                    add.params = list(color = "red",
                                                      fill = "lightgray"),
                                    palette = "jco",     # journal color palette. see ?ggpar
                                    xlab = x.label,      # x-axis label
                                    ylab = "SMR",        # y-axis label
                                    title = plot.title)+    
                  ggpubr::stat_cor(label.y = 5)+         # Add correlation coefficient
                  ggpubr::font("title", size = 10, face = "bold")+
                  ggpubr::font("xlab", size = 10)+
                  ggpubr::font("ylab", size = 10)
 
  return(scatter.plot)
  
}

# Invoking function 
# 2017 scatter plots ------------------------------------------------------------

scatter_1km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_1km",
                                       x.label = "Percentage of catchment population \nliving in 1km buffer",
                                       plot.title = "2017")

scatter_2km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_2km",
                                        x.label = "Percentage of catchment population \nliving in 2km buffer",
                                        plot.title = "2017")

scatter_3km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_3km",
                                        x.label = "Percentage of catchment population \nliving in 3km buffer",
                                        plot.title = "2017")

# 2018 scatter plots -----------------------------------------------------------
scatter_1km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_1km",
                                        x.label = "Percentage of catchment population \nliving in 1km buffer",
                                        plot.title = "2018")

scatter_2km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_2km",
                                        x.label = "Percentage of catchment population \nliving in 2km buffer",
                                        plot.title = "2018")

scatter_3km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_3km",
                                        x.label = "Percentage of catchment population \nliving in 3km buffer",
                                        plot.title = "2018")

# 2019 scatter plots -----------------------------------------------------------
scatter_1km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_1km",
                                        x.label = "Percentage of catchment population \nliving in 1km buffer",
                                        plot.title = "2019")

scatter_2km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_2km",
                                        x.label = "Percentage of catchment population \nliving in 2km buffer",
                                        plot.title = "2019")

scatter_3km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_3km",
                                        x.label = "Percentage of catchment population \nliving in 3km buffer",
                                        plot.title = "2019")

# 2020 scatter plots -----------------------------------------------------------
scatter_1km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_1km",
                                        x.label = "Percentage of catchment population \nliving in 1km buffer",
                                        plot.title = "2020")

scatter_2km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_2km",
                                        x.label = "Percentage of catchment population \nliving in 2km buffer",
                                        plot.title = "2020")

scatter_3km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_3km",
                                        x.label = "Percentage of catchment population \nliving in 3km buffer",
                                        plot.title = "2020")

# Arrange the plots ------------------------------------------------------------
ggpubr::ggarrange(scatter_1km_2017, scatter_2km_2017, scatter_3km_2017,
                  scatter_1km_2018, scatter_2km_2018, scatter_3km_2019,
                  scatter_1km_2019, scatter_2km_2019, scatter_3km_2019, 
                  scatter_1km_2020, scatter_2km_2020, scatter_3km_2020,
                  ncol = 3, nrow = 4)
Fig 10. Relationship between standardised morbidity ratio and living near waterbodies

Fig 10. Relationship between standardised morbidity ratio and living near waterbodies

Model fitting

# Combine data for model fitting -----------------------------------------------
model_data_2017 <- merge(expected_malaria_2017, smr_pop_2017, by = "rowID", all = TRUE) %>% 
  dplyr::select(-Names.y) %>% 
  dplyr::rename(Names = Names.x)

model_data_2018 <-  merge(expected_malaria_2018, smr_pop_2018, by = "rowID", all = TRUE) %>% 
  dplyr::select(-Names.y) %>% 
  dplyr::rename(Names = Names.x)

model_data_2019 <-  merge(expected_malaria_2019, smr_pop_2019, by = "rowID", all = TRUE) %>% 
  dplyr::select(-Names.y) %>% 
  dplyr::rename(Names = Names.x)

model_data_2020 <-  merge(expected_malaria_2020, smr_pop_2020, by = "rowID", all = TRUE) %>% 
  dplyr::select(-Names.y) %>% 
  dplyr::rename(Names = Names.x)

# Fit generalised linear model -------------------------------------------------
# Defining model parameters:
# response variable: observed_2017, observed_2018, observed_2019, observed_2020 are
#                    recorded dry season malaria cases in that year
# risk factor: prop_pop_1km,  prop_pop_2km,  prop_pop_3km are the percentage of people living
#             within 1km, 2km and 3km buffers of water bodies, respectively.
# offset: expected_year is the number of malaria cases we would expect if the malaria rate
#          was equal in all the catchment areas

# 2017 -------------------------------------------------------------------------
model_1km_2017 <- glm(observed_2017~1+prop_pop_1km+offset(log(expected_2017)),
                      data = model_data_2017, family = 'poisson')

summary(model_1km_2017)
## 
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_1km + offset(log(expected_2017)), 
##     family = "poisson", data = model_data_2017)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -50.150  -19.463   -1.324   20.443  118.758  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2223780  0.0059167  -37.58   <2e-16 ***
## prop_pop_1km  0.0247164  0.0004965   49.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33562  on 23  degrees of freedom
## Residual deviance: 31130  on 22  degrees of freedom
##   (3 observations deleted due to missingness)
## AIC: 31372
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_1km_2017)
  observed_2017
Predictors Incidence Rate Ratios CI p
(Intercept) 0.80 0.79 – 0.81 <0.001
prop_pop_1km 1.03 1.02 – 1.03 <0.001
Observations 24
R2 Nagelkerke 1.000
model_2km_2017 <- glm(observed_2017~1+prop_pop_2km+offset(log(expected_2017)),
                      data = model_data_2017, family = 'poisson')

summary(model_2km_2017)
## 
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_2km + offset(log(expected_2017)), 
##     family = "poisson", data = model_data_2017)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -45.875  -19.972    1.747   19.501  118.324  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3126953  0.0062855  -49.75   <2e-16 ***
## prop_pop_2km  0.0134272  0.0002117   63.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33564  on 24  degrees of freedom
## Residual deviance: 29695  on 23  degrees of freedom
##   (2 observations deleted due to missingness)
## AIC: 29947
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_2km_2017)
  observed_2017
Predictors Incidence Rate Ratios CI p
(Intercept) 0.73 0.72 – 0.74 <0.001
prop_pop_2km 1.01 1.01 – 1.01 <0.001
Observations 25
R2 Nagelkerke 1.000
model_3km_2017 <- glm(observed_2017~1+prop_pop_3km+offset(log(expected_2017)),
                      data = model_data_2017, family = 'poisson')

summary(model_3km_2017)
## 
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_3km + offset(log(expected_2017)), 
##     family = "poisson", data = model_data_2017)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -47.915  -21.148    1.191   19.098  121.065  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3936999  0.0075710  -52.00   <2e-16 ***
## prop_pop_3km  0.0100819  0.0001629   61.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33564  on 24  degrees of freedom
## Residual deviance: 29726  on 23  degrees of freedom
##   (2 observations deleted due to missingness)
## AIC: 29977
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_3km_2017)
  observed_2017
Predictors Incidence Rate Ratios CI p
(Intercept) 0.67 0.66 – 0.68 <0.001
prop_pop_3km 1.01 1.01 – 1.01 <0.001
Observations 25
R2 Nagelkerke 1.000
# 2018 -------------------------------------------------------------------------
model_1km_2018 <- glm(observed_2018~1+prop_pop_1km+offset(log(expected_2018)),
                      data = model_data_2018, family = 'poisson')

summary(model_1km_2018)
## 
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_1km + offset(log(expected_2018)), 
##     family = "poisson", data = model_data_2018)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -55.722  -12.556    3.871   13.940  137.674  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3987643  0.0060731  -65.66   <2e-16 ***
## prop_pop_1km  0.0323262  0.0004041   79.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 42339  on 25  degrees of freedom
## Residual deviance: 35921  on 24  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 36187
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_1km_2018)
  observed_2018
Predictors Incidence Rate Ratios CI p
(Intercept) 0.67 0.66 – 0.68 <0.001
prop_pop_1km 1.03 1.03 – 1.03 <0.001
Observations 26
R2 Nagelkerke 1.000
model_2km_2018 <- glm(observed_2018~1+prop_pop_2km+offset(log(expected_2018)),
                      data = model_data_2018, family = 'poisson')

summary(model_2km_2018)
## 
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_2km + offset(log(expected_2018)), 
##     family = "poisson", data = model_data_2018)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -58.860  -11.691    4.018   20.409  114.785  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.5633136  0.0067344  -83.65   <2e-16 ***
## prop_pop_2km  0.0167495  0.0001679   99.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 42339  on 25  degrees of freedom
## Residual deviance: 32218  on 24  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 32484
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_2km_2018)
  observed_2018
Predictors Incidence Rate Ratios CI p
(Intercept) 0.57 0.56 – 0.58 <0.001
prop_pop_2km 1.02 1.02 – 1.02 <0.001
Observations 26
R2 Nagelkerke 1.000
model_3km_2018 <- glm(observed_2018~1+prop_pop_3km+offset(log(expected_2018)),
                      data = model_data_2018, family = 'poisson')

summary(model_3km_2018)
## 
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_3km + offset(log(expected_2018)), 
##     family = "poisson", data = model_data_2018)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -52.920  -14.159    5.671   22.270  106.698  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.8634330  0.0089250  -96.74   <2e-16 ***
## prop_pop_3km  0.0165420  0.0001522  108.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 42339  on 25  degrees of freedom
## Residual deviance: 29688  on 24  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 29955
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_3km_2018)
  observed_2018
Predictors Incidence Rate Ratios CI p
(Intercept) 0.42 0.41 – 0.43 <0.001
prop_pop_3km 1.02 1.02 – 1.02 <0.001
Observations 26
R2 Nagelkerke 1.000
# 2019 -------------------------------------------------------------------------
model_1km_2019 <- glm(observed_2019~1+prop_pop_1km+offset(log(expected_2019)),
                      data = model_data_2019, family = 'poisson')

summary(model_1km_2019)
## 
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_1km + offset(log(expected_2019)), 
##     family = "poisson", data = model_data_2019)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -66.815  -18.307    1.353   23.118  102.556  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3353561  0.0077298  -43.38   <2e-16 ***
## prop_pop_1km  0.0180337  0.0003618   49.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33446  on 26  degrees of freedom
## Residual deviance: 30926  on 25  degrees of freedom
## AIC: 31194
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_1km_2019)
  observed_2019
Predictors Incidence Rate Ratios CI p
(Intercept) 0.72 0.70 – 0.73 <0.001
prop_pop_1km 1.02 1.02 – 1.02 <0.001
Observations 27
R2 Nagelkerke 1.000
model_2km_2019 <- glm(observed_2019~1+prop_pop_2km+offset(log(expected_2019)),
                      data = model_data_2019, family = 'poisson')

summary(model_2km_2019)
## 
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_2km + offset(log(expected_2019)), 
##     family = "poisson", data = model_data_2019)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -62.833  -21.458    2.378   24.375  101.995  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3051901  0.0088349  -34.54   <2e-16 ***
## prop_pop_2km  0.0066732  0.0001748   38.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33446  on 26  degrees of freedom
## Residual deviance: 31955  on 25  degrees of freedom
## AIC: 32223
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_2km_2019)
  observed_2019
Predictors Incidence Rate Ratios CI p
(Intercept) 0.74 0.72 – 0.75 <0.001
prop_pop_2km 1.01 1.01 – 1.01 <0.001
Observations 27
R2 Nagelkerke 1.000
model_3km_2019 <- glm(observed_2019~1+prop_pop_3km+offset(log(expected_2019)),
                      data = model_data_2019, family = 'poisson')

summary(model_3km_2019)
## 
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_3km + offset(log(expected_2019)), 
##     family = "poisson", data = model_data_2019)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -60.344  -22.596    4.515   23.613  102.074  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3214559  0.0105946  -30.34   <2e-16 ***
## prop_pop_3km  0.0046685  0.0001437   32.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33446  on 26  degrees of freedom
## Residual deviance: 32357  on 25  degrees of freedom
## AIC: 32624
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_3km_2019)
  observed_2019
Predictors Incidence Rate Ratios CI p
(Intercept) 0.73 0.71 – 0.74 <0.001
prop_pop_3km 1.00 1.00 – 1.00 <0.001
Observations 27
R2 Nagelkerke 1.000
# 2020 -------------------------------------------------------------------------
model_1km_2020 <- glm(observed_2020~1+prop_pop_1km+offset(log(expected_2020)),
                      data = model_data_2020, family = 'poisson')

summary(model_1km_2020)
## 
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_1km + offset(log(expected_2020)), 
##     family = "poisson", data = model_data_2020)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -83.14  -23.49    9.60   24.68   89.43  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.0236328  0.0060477  -3.908 9.32e-05 ***
## prop_pop_1km  0.0024199  0.0004792   5.050 4.41e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 36512  on 22  degrees of freedom
## Residual deviance: 36486  on 21  degrees of freedom
##   (4 observations deleted due to missingness)
## AIC: 36719
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_1km_2020)
  observed_2020
Predictors Incidence Rate Ratios CI p
(Intercept) 0.98 0.97 – 0.99 <0.001
prop_pop_1km 1.00 1.00 – 1.00 <0.001
Observations 23
R2 Nagelkerke 0.670
model_2km_2020 <- glm(observed_2020~1+prop_pop_2km+offset(log(expected_2020)),
                      data = model_data_2020, family = 'poisson')

summary(model_2km_2020)
## 
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_2km + offset(log(expected_2020)), 
##     family = "poisson", data = model_data_2020)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -79.646  -30.318    7.777   26.260   87.868  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.1690115  0.0065167  -25.93   <2e-16 ***
## prop_pop_2km  0.0056305  0.0002104   26.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 37822  on 23  degrees of freedom
## Residual deviance: 37108  on 22  degrees of freedom
##   (3 observations deleted due to missingness)
## AIC: 37350
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_2km_2020)
  observed_2020
Predictors Incidence Rate Ratios CI p
(Intercept) 0.84 0.83 – 0.86 <0.001
prop_pop_2km 1.01 1.01 – 1.01 <0.001
Observations 24
R2 Nagelkerke 1.000
model_3km_2020 <- glm(observed_2020~1+prop_pop_3km+offset(log(expected_2020)),
                      data = model_data_2020, family = 'poisson')

summary(model_3km_2020)
## 
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_3km + offset(log(expected_2020)), 
##     family = "poisson", data = model_data_2020)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -78.13  -29.63   10.56   25.70   86.49  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2336361  0.0075204  -31.07   <2e-16 ***
## prop_pop_3km  0.0048908  0.0001526   32.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 37822  on 23  degrees of freedom
## Residual deviance: 36775  on 22  degrees of freedom
##   (3 observations deleted due to missingness)
## AIC: 37018
## 
## Number of Fisher Scoring iterations: 5
sjPlot::tab_model(model_3km_2020)
  observed_2020
Predictors Incidence Rate Ratios CI p
(Intercept) 0.79 0.78 – 0.80 <0.001
prop_pop_3km 1.00 1.00 – 1.01 <0.001
Observations 24
R2 Nagelkerke 1.000

Check how well the fitted values line up with the observations

The fitted values appear to line up particularly well with the observed data, suggesting that prop_pop_* i.e., proportion of catchment population living near water bodies can help us understand malaria risk in the catchment areas.

# Helper function to create scatter plots to see how well 
# fitted values line up with observed malaria cases
plot.fitted.values <- function(fitted.values.df, model.df, title){
  
  # Remove missing values from model data since 
  # model fitting deletes missing observations
  model.df.complete <- model.df %>% 
    tidyr::drop_na() %>% 
    dplyr::rename_at(vars(starts_with("observed_")), ~ str_c("observed"))

  # Plot fitted versus observed values
  scatter.plot <- ggplot2::ggplot()+ 
                  ggplot2::geom_point(aes(fitted.values.df$fitted.values, 
                                      model.df.complete$observed))+
                  ggplot2::theme_classic()+
                  ggplot2::labs(y = "Observed values", 
                                x = "Fitted values",
                                title = title)
  return(scatter.plot)
}

# Invoking function 
# 2017 -------------------------------------------------------------------------
fitted_1km_2017 <- plot.fitted.values(model_1km_2017, model_data_2017, "2017: 1km model")

fitted_2km_2017 <- plot.fitted.values(model_2km_2017, model_data_2017, "2017: 2km model")

fitted_3km_2017 <- plot.fitted.values(model_3km_2017, model_data_2017, "2017: 3km model")

# 2018 -------------------------------------------------------------------------
fitted_1km_2018 <- plot.fitted.values(model_1km_2018, model_data_2018, "2018: 1km model")

fitted_2km_2018 <- plot.fitted.values(model_2km_2018, model_data_2018, "2018: 2km model")

fitted_3km_2018 <- plot.fitted.values(model_3km_2018, model_data_2018, "2018: 3km model")

# 2019 -------------------------------------------------------------------------
fitted_1km_2019 <- plot.fitted.values(model_1km_2019, model_data_2019, "2019: 1km model")

fitted_2km_2019 <- plot.fitted.values(model_2km_2019, model_data_2019, "2019: 2km model")

fitted_3km_2019 <- plot.fitted.values(model_3km_2019, model_data_2020, "2019: 3km model")

# 2020 -------------------------------------------------------------------------
fitted_1km_2020 <- plot.fitted.values(model_1km_2020, model_data_2020, "2020: 1km model")

fitted_2km_2020 <- plot.fitted.values(model_2km_2020, model_data_2020, "2020: 2km model")

fitted_3km_2020 <- plot.fitted.values(model_3km_2020, model_data_2020, "2020: 3km model")

# Layout scatter plots ---------------------------------------------------------
cowplot::plot_grid(fitted_1km_2017, fitted_1km_2018, fitted_2km_2018, fitted_3km_2018,
                   fitted_1km_2019, fitted_2km_2019, fitted_1km_2020, ncol = 3, nrow = 4)
Fig 10. How well the percentage of catchment population living around water bodies explain observed malaria incidence

Fig 10. How well the percentage of catchment population living around water bodies explain observed malaria incidence

Test for residual spatial autocorrelation

Here we’ll use adjacency as criterion

# Find adjacent polygons,